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Significant evidence for a gravitational-wave background was reported by several pulsar- 
timing-array collaborations. By assuming that this signal is interpreted by the scalar-induced 
gravitational waves, we study physical implications of the observed signal for the nature of 
primordial curvature perturbations and primordial black holes. In particular, we explore the 
effects of primordial non-Gaussianity on the inferences of model parameters, and obtain the 
parameter region allowed by the observed signal, i.e., the primordial scalar spectral amplitude 
Ag ~ 107? — 1, the primordial non-Gaussian parameter —10 < fn < 10, and the mass of 
primordial black holes mppn ~ 1078 —0.1Mo. We find that the non-Gaussianity suppressing 
the abundance of primordial black holes is preferred by the observed signal. We show that 
the anisotropies of scalar-induced gravitational waves are a powerful probe for measurements 
of the non-Gaussian parameter fnL, and conduct a complete analysis of the angular power 
spectrum in the nano-Hertz band. We expect that the Square Kilometre Array project has 


potentials to measure such anisotropies. 


I. INTRODUCTION 


Recently, several observational collaborations of pulsar timing array (PTA) reported signifi- 
cant evidence for an excess signal with Hellings-Downs (HD) correlations [1—4], indicating the 
gravitational-wave origin of this signal. The strain amplitude of such a gravitational-wave back- 
ground (GWB) was found to be of order 10715 at the pivot frequency 1 yr~!. The inferred GWB 
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spectrum was found to be consistent with the astrophysical origin due to inspiraling supper-massive 
black hole (SMBH) binaries [5]. However, the current datasets could not exclude possibilities of 
cosmological origins (and other exotic astrophysical sources), which were studied by the collabora- 
tions in several accompany papers [6, 7]. In particular, many models of cosmological origins have 
been shown to provide even better fits to the observed signal than the SMBH-binary interpretation. 
If confirmed in future, they may point to evidence for new physics. 

In this work, we will focus on the cosmological interpretation to the observed signal, i.e., the 
scalar-induced gravitational waves (SIGWs) [8-13]. This possibility was previously considered to 
account for the NANOGrav 12.5-year dataset [14] by the authors of Refs. [15-24]. Recently, it 
was revisited by the collaborations in Refs. [6, 7], but only the Gaussian primordial scalar (or 
equivalently, comoving curvature) perturbations were considered in the above studies. However, 
it has been shown that the primordial non-Gaussianity contributes a lot to the energy density of 
SIGWs [25-33], indicating significant modifications to the energy-density fraction spectrum that 
is essential for the data analysis of PTA observations. By interpreting the observed signal to 
have a SIGW origin, therefore, we will study implications of the PTA datasets for the nature of 
primordial scalar perturbations, including the power spectrum and the local-type primordial scalar 
non-Gaussianity. 

We will also study physical implications for scenarios of primordial black holes (PBHs). The 
formation process of PBHs was accompanied by the inevitable production of SIGWs. In fact, the 
enhanced curvature perturbations not only produced PBHs due to gravitational collapse in the 
early universe [34], but also induced a GWB via nonlinear mode-couplings. Therefore, we can 
explore the PBH scenario by using SIGWs [35-38]. Related works analyzing realistic datasets can 
be found in Refs. [6, 7, 20-22, 38, 39]. Some strong influences of primordial non-Gaussianity on 
the mass function of PBHs were also studied in literature [24, 40-52]. In this work, by taking 
into account the effects of primordial non-Gaussianity, we will recast the constraints on primordial 
perturbations into the constraints on the mass function of PBHs '. 

To further explore the primordial non-Gaussiantiy, we will study the anisotropies in SIGWs in 
the PTA band and provide a complete analysis for the angular power spectrum in the same band. 
It was known that the energy-density spectrum of isotropic SIGWs has significant degeneracies of 


1 Note added: At the time of writing the current paper, there appeared a related paper [53] that analyzed the 
posteriors of NANOGrav 15-year (NG15) data. The authors claimed that the Gaussian scenarios for SIGWs are 
in tension with the current PTA data at 20 confidence level, but the non-Gaussian scenarios that suppress the 
abundance of PBHs can alleviate such a tension. Due to large uncertainties in formation scenarios of PBHs (e.g., 
see reviews in Ref. [54]), one of the leading aims for our current work is a comprehensive study on the primordial 
non-Gaussianity. Full Bayesian analysis taking into the effects of primordial non-Gaussianity on SIGWs will be 


left to our future works. 


model parameters [33]. Conducting a complete analysis of angular power spectrum, we showed that 
the angular power spectrum is useful for determining fn, because the above degeneracies would 
be broken if the anisotropic SIGWs are considered [33]. Earlier related works can be found in 
Refs. [54-62]. In this work, in the PTA band, we will further study the angular power spectrum for 
the anisotropies in SIGWs, which may be not only useful for determining fy , but also important 
for discriminating different GWB sources in realistic data analysis. 

The remaining context of this paper is arranged as follows. In Section lI, we will provide a brief 
summary of the isotropic SIGWs. In Section III, we will show implications of the current datasets 
for the power spectrum of primordial curvature perturbations and then for the mass function 
of PBHs. In Section IV, we will study the anisotropies in SIGWs and show the angular power 


spectrum in PTA band. In Section V, we make concluding remarks. 


IL. ENERGY-DENSITY FRACTION SPECTRUM OF SCALAR-INDUCED 
GRAVITATIONAL WAVES 


In this section, we show a brief but self-consistent summary of the main results of SIGW theory. 

The energy-density fraction spectrum of the isotropic GWB is Oge(7, q) = Pgew(N, @)/ Perit (7) 
[63], where q is the wavenumber of gravitational waves (GWs), pcrit is the critical energy density 
of the universe at the conformal time 7, and the overbar denotes physical quantities of background 
level. The above definition indicates that f pew(7,q)dlnq is the total energy-density fraction of 
GWB [63]. The spectrum is formally expressed as Pgw(n, q) ~ (hijthiza), Where hij(n, q) denotes 
the strain with wavevector q in Fourier space and the angle brackets define an ensemble average. 
For SIGWs on subhorizon scales, we have hij ~ ¢? and then Qgw(7,q) ~ (C4) [8, 9], where ¢(q) 
denotes the primordial curvature perturbations in the early universe. In the case of primordial 
Gaussianity, the semi-analytic formula for Qgw(n, q) was shown in Refs. [12, 13] and earlier related 
works can be found in Refs. [8, 9]. However, in the case of primordial non-Gaussianity, there is not 
such a semi-analytic formula. Related studies can be found in recent literature [25-33, 64]. In this 
work, we will follow the conventions of our previous work [33]. 

To derive the contributions of primordial non-Gaussianity to the energy density, we express the 


primordial curvature perturbations ¢ in terms of their Gaussian components Cg, i.e., [65] 
3 


3 d°’k 
la) = Gla) + Fx | geala- k) , (1) 
where fni is the non-linear parameter characterizing the local-type primordial non-Gaussianity. 


We introduce a new quantity of Fy, = 3fnL/5 that will simplifies the analytic formulae in the 
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FIG. 1. Unscaled (or equivalently, Ag = 1 and Fy, = 1) contributions to the energy-density fraction 
spectrum of isotropic SIGWs. We produce this figure by using the original data of Ref. [33]. 


following. Note that the validation of perturbation theory requires Ag FR, < 1, where Ag will be 


defined in the following. We define the dimensionless power spectrum of ¢, as follows 


T? 
lla la’) = O(a + a) ARM) , (2) 


where AQ) is assumed to be a normal function with respect to ln q in this work, i.e., [21, 31, 66, 67] 
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Here, Ag stands for the spectral amplitude at the spectral peak wavenumber q,, and ø denotes the 
standard deviation characterizing the width of spectrum. The wavenumber q can be straightfor- 
wardly recast into the frequency v, namely, q = 27v and q = 271. 

By conducting a tedious but straightforward derivation process, we can decompose Oi ~ (C4) 
into three components depending on the power of ft, based on the Wick’s theorem. However, 
the complete derivations have been simplified by following an approach of Feynman-like diagrams 


[25, 28, 30-33]. We summarize only the final results as follows 


Qew(n, a) = On, a) + OQ (ng) + A(n, a), (4) 


where we show the analytic expressions of Qe x AZ (As Fl)” explicitly in Appendix A. They 


have been computed with the vegas package [68], and the numerical results are reproduced in 
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FIG. 2. Energy-density fraction spectrum of isotropic SIGWs. The NANOGrav 15-year data points [7] are 


shown for comparison. Here, we fix o = 1. 


Fig. 1. In particular, a is exactly the result of energy-density fraction spectrum corresponding 


to the case of primordial Gaussianity, while AW and a completely describe the contributions of 
local-type primordial non-Gaussianity. 


The energy-density fraction spectrum of SIGWs at the current conformal time no is given by 


4/3 
cen E T (2) Gece mcr (5) 


where Qraaoh? = 4.2 x 10~° is the physical energy-density fraction of radiations in the present 
universe [69], T (and Teq) labels the cosmic temperatures at the emission time (and the epoch of 


matter-radiation equality), and v denotes the gravitational-wave frequency, i.e., [21] 


1/2 —1/3 
nHz GeV 106.75 106.75 


Here, the effective relativistic degrees of the universe, i.e., gx p and gx s, are tabulated functions of 


T, as shown in Ref. [70]. To illustrate the SIGW-interpretation of the current datasets, we depict 


Qew,o(v) as a function of v in Fig. 2, by choosing a particular set of model parameters. 


II. IMPLICATIONS OF PTA DATASETS FOR PRIMORDIAL CURVATURE 
PERTURBATIONS AND PRIMORDIAL BLACK HOLES 


In this section, we study possible constraints on the parameter space of primordial power spec- 
trum and PBHs from the NG15 data. Constraints from other PTA datasets can also be obtained 


following the same approach, but disregarded in this work. 
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FIG. 3. Contours of the 68% probability regions for the spectral amplitude Ag and peak frequency vs 
defined in Eq. (3). Different values of |FyL| have been considered, but ø is fixed to unity. 


A. Primordial curvature perturbations 


Based on the principle of PTA observations, the energy density of a given GWB, denoted with 


Qew(v) here, is related with the timing residual power spectral density S(v), i.e., [7] 


Srt’ v i 
Qge(¥) = = (=) S(v) , (7) 


where Vyr is a pivot frequency related to a duration time of one year. For the realistic data analysis, 


S(v) could be assumed to be power-law, i.e., [7] 


ar 
es (=) yi, (8) 


~ 1072 Vyr 


where the amplitude A and index y have been constrained by the NG15 data. For example, based 
on Fig. 1 of Ref. [7], the green contours in log}, A — y plane stand for 68% (inner) and 95% (outer) 
probability regions. 
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FIG. 4. Validation of perturbation theory versus the inferred parameter region (68% confidence level) from 
the NG15 data. The dotted lines denote Fyt = —(4¢.)~! and the solid curves denote mppn = 107? Mo and 
fobn = 1 in the case of Çe = 0.7 (purple color) and Çe = 1.2 (rose color). 


Though a full Bayesian analysis is of necessity in the future, we can immediately get some useful 
insights by recasting the contours in log}, A — y plane into the contours in Ag — v, plane, once 
o and |FN | take definitive values. Fig. 3 displays such contours for a series of values of |FNi], 
denoted with colored dots. Hereafter, we fix the value of o, namely, o = 1, but a generalization is 
straightforward. In fact, such a choice is consistent with the best-fit o from the NG15 data, which 
resulted in a 68% credible interval ø € [0.51, 2.07] [7]. To conduct such recasting manipulations, for 
a given set of A and y, we should vary Ag and v, to make New olv) following the same power-law 
as Qew(v). Here, we have fixed the pivot scale to be Vyr for the model parameters Ag and v4. The 
above approach was already adopted in literature [17]. 

Based on Fig. 3, we find significant effects of non-Gaussian parameter |Fyz| on the inferences 
of other parameters Ag and ną, and vice versa. In particular, the contours are shifted to lower-Ag 
regimes with the increase of | FNL], since the energy density of SIGWs would be overproduced in the 
presence of primordial non-Gaussianity. This is one of the important results of this work. Though 
our research method can be straightforwardly generalized to other values of ø, it is essential to 
conduct a full Bayesian analysis that is left to our future works. Note that the above results are 
independent of the non-Gaussian sign. However, the non-Gaussian sign would significantly affect 


the mass function of PBHs that will be studied in the following. 


We further pin down the parameter region of Ag and FNL, by considering the validation of 
perturbation theory that requires Ask, < 1. In Fig. 4, the blue shaded area corresponds to such 
a requirement, while the orange shaded area corresponds to the inferred parameter region (68% 
confidence level) from the NG15 data. Only the overlap area is simultaneously allowed by theory 
and observations. In this sense, the PTA observations have already been a powerful probe of the 


early universe. 


B. Primordial black holes 


We further recast the constraints on the primordial curvature perturbations into the constraints 
on PBHs. Due to large uncertainties for the formation scenarios of PBHs (e.g., see reviews in 
Ref. [54]), we consider a simplified scenario used by Ref. [48] to demonstrate the importance of 


primordial non-Gaussianity. The initial mass function of PBHs is given by 


g=[  P@ac= f : ep ( i) da, (9) 
oe C(eg)>Ce W210 g 20? 


where P(Ç) is a probability distribution function (PDF) of primordial curvature perturbations, oy 


is a standard variance for a PDF of Gaussian component Çg, and Çe denotes the critical fluctuation. 
Considering the power spectrum in Eq. (3), we obtain 07 = (¢?) = f A2(q)dIng = As. In addition, 
it is known that Çe ~ O(1), which was shown to be 0.7 and 1.2 in Ref. [71]. We will consider both 
of them in the following. 

To calculate Eq. (9), we separate Fyi into two regimes, i.e., Fy, > 0 and Fyni < 0. As the first 


step, we solve the equation ¢(¢,) = Çe to get the following relation 


a= —-14+VJ71+4+ 4FNLG (10) 
gx 2FNL , 


By substituting it into Eq. (9), we have an expression of 3 for Fyni > 0, i.e., 


= s-( fË af I (Cg)dlg = serfe (4) + zefe( E) , (11) 


where erfc(x) is the complementary error function. Similarly, for —(4¢.)~' < Fyi < 0, we have 


p= ie P(Cy)dCy = serfe a serfe cS) (12) 


However, for Fy, < —(4¢-)~!, no PBHs were formed in the early universe, because the curvature 


perturbations are expected to never exceed the critical fluctuation. As a reasonable candidate of 


cold dark matter, the abundance of PBHs is given as [72] 


—1/4 —1/2 
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where Mpph Stands for the mass of PBHs in units of Mo, and Tẹ denotes the cosmic temperature 
at the formation time. Roughly speaking, mppn can be related with the horizon mass my and then 
the frequency v. To be specific, we have [17] 


to E (14 
Mo ~ 0.31Mo T 5.0nHz ` 


Therefore, we infer the mass of PBHs to be ~ O(107? — 107t)Mọo, based on Fig. 3. However, 
the inferred abundance would be larger than unity, indicating that the PBH scenario is in tension 
with the NG15 data. To demonstrate this result more clearly, in Fig. 4, we depict two solid curves 
corresponding to Mpbh = 107? Mo and fobh = 1 in the case of Çe = 0.7 (purple curve) and ¢ = 1.2 
(rose curve). We also denote the critical value Fyr = —(4¢.)~! with vertical dotted lines. 

Based on Fig. 4, when we interpret the observed signal to have a SIGW origin, we find overpro- 
duction of PBHs, since the inferred value of Ag is typically one order of magnitude larger than the 
value of Ag producing fpbh = 1. Even if we take into account the effects of positive non-Gaussianity, 
such overproduction can not be effectively alleviated. In contrast, the negative non-Gaussianity 
can alleviate it, particularly when we consider a sizable negative non-Gaussian parameter, i.e., 
Fy, < —(4¢.)~!, that forbids any formation of PBHs. Due to large uncertainties, it is challenging 
to exclude the PBH scenario by analyzing the current datasets. However, a more-detailed analysis 
such as the Bayesian analysis is still important, but beyond the scope of the current paper. 

Therefore, it seems essential to measure the primordial non-Gaussianity, at least determine 
the sign of Fy in order to judge the PBH scenario. However, due to the sign degeneracy of 
Fy, it is impossible to determine the non-Gaussian sign via measurements of the energy-density 
fraction spectrum of SIGWs. In the next section, we will propose that the anisotropic SIGWs 
have potentials to break such a sign degeneracy as well as other degeneracies of model parameters, 


providing some possibilities for judgements of the PBH scenario in future. 


IV. ANGULAR POWER SPECTRUM FOR ANISOTROPIES IN SCALAR-INDUCED 
GRAVITATIONAL WAVES 


In this section, we study the anisotropies in SIGWs as well as the angular power spectrum in 
the PTA band, following the conventions of our previous paper [33]. 

The anisotropies in SIGWs arise from long-wavelength modulations of the energy density pro- 
duced by short-wavelength modes. Based on Section II, we know that SIGW were produced at 


extremely high redshifts, corresponding to extremely small horizons. Due to limitation in the 
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angular resolution of gravitational-wave detectors, the signal along a line-of-sight stands for an 
ensemble average of the energy densities over a quantity of such horizons. In this sense, any two 
signals would be identical. However, the energy density of SIGWs that were produced by short- 
wavelength modes can be spatially redistributed by long-wavelength modes, if there are couplings 
between the short- and long-wavelength modes. Therefore, the primordial non-Gaussianity of local 
type could have contributions to such couplings, as shown in the following. 

Analogue to the temperature fluctuations of relic photons [73], the initial inhomogeneities 
in SIGWs at spatial location x are characterized with the density contrast, i.e., dsw(7,x,q) = 
AtWwew(n, X;Q)/Qew(,q) — 1. Here, the energy-density full spectrum wew(n,x,q) is defined by 
the energy density pgw(7,x) = Perit J d°q Wew(7,x,q)/g?. This definition implies that we have 
Wew(7,X,q) ~ (¢*)x, where the subscript x stands for an ensemble average within the horizon 
enclosing x [33, 55]. We further separate ¢, into short-wavelength gs and long-wavelength Cyr, 
i.e., Cg = Cgs + gr [74]. Finally, we obtain dgw(7,x,q) ~ Cg (Cgs¢2)x at linear order of Çgz. Here, 
we introduce Çs to denote the part of ¢, which is composed of only gs. Higher orders of gz are 
negligible due to power spectrum AZ ~ 107° for Cgz, [69]. By using the Feynman-like rules and 
diagrams, we have obtained an explicit expression for dgw(7,x,q) as follows [33] 


3 fNL a dk 
5 Qgw(7, 9) (2m)3/2 


where a new quantity is introduced for simplifying our computation, i.e., 


dew(, x, q) = e “Cot (k ) ’ (15) 


Qne(n, a) = BAO (n, q) + ZARN, a) - (16) 


The initial inhomogeneities here correspond to the intrinsic temperature fluctuations of cosmic 
microwave background (CMB) photons on the last-scattering surface. 

The “observed” density contrast sw, o(q) can be analytically estimated following the line-of-sight 
approach [75-77]. It is contributed by both of the initial inhomogeneities and some propagation 
effects, i.e, [33] 


Sgw,0(G) = ôgw(1, x, a) + [4 — New,o(¥)] (7; %) ; (17) 


where we consider only the Sachs-Wolfe (SW) effect [78], characterized by the Bardeen’s potential 


at large scales, i.e., 


ik: 
1), x =. JE 27) ae e *CgL(k) , (18) 
and Ngw,o(q) is the index of the energy-density fraction spectrum defined in Eq. (5), i.e., 


iais Oln Oggi) _ ln a 
ae Olnv 7 Olng q=2rv 


(19) 
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In the following, we adopt the assumption of statistical isotropy for the density contrasts on large 
scales. Analogue to the study of CMB (e.g., see Ref. [79]), the inhomogeneities in SIGWs can be 
recast into the anisotropies in GWB. 

The reduced angular power spectrum is usually used for characterizing the statistics of the 


anisotropies in SIGWs. It is defined as the following two-point correlator 
(Sgw,0,6m (27) 5w0,e%m’ (27) = ðw ômm Ce(V) , (20) 


where dgw,o(q) is expanded in terms of spherical harmonics, i.e., dgw,o(@) = oom Sew,0,m(G) Yom (0). 
Roughly, we get Cy ~ Ori x (CgLCgL) ~ AZ. A complete analysis employing the Feynman-like 
rules and diagrams has been conducted in our previous work [33]. We summarize only the final 


results as follows 


2 


T 2 
2e tasr (21) 


= Qng(N, 2rv) 
~ -25£(€+ 1) 


fNu= 


Cel) Qew(n, 27) 


which can be recast into the angular power spectrum, i.e., 


2 
C= (Pea) Cr(v) . (22) 


We make an analogue to the anisotropies in CMB, for which the rms temperature is roughly de- 
termined by [¢(¢+1)C£M®/(27)|!/?. For the anisotropies in SIGWs, the rms energy density is 
roughly determined by [é(€+ 1)Ce(v)/(27)]'/?, which is the variance of the energy-density fluctu- 
ations. Note that the rms energy density is constant with respect to the multipoles £, but depends 
on the gravitational-wave frequency bands. 

In Fig. 5, we depict the rms energy density with respect to the gravitational-wave frequency. 
For comparison, we also show the energy-density fraction spectrum in the same figure. Roughly 
speaking, we find Ver ~ O(10~*), depending on the model parameters. Note that the angular 
power spectrum can break the degeneracies of model parameters. For example, based on Fig. 5, 
we find coincidence for the energy-density fraction spectra of three different sets of parameters. 
However, the angular power spectrum breaks such a coincidence, particularly, the sign degeneracy 
of Fyy. This result indicates that the primordial non-Gaussianity could be determined by mea- 
surements of the anisotropies in SIGWs in principle [33]. Recently, an upper limit on the reduced 
angular power spectrum was inferred to be Cp < 20% from the NG15 data [81]. It is not precise 
enough to test the theoretical predictions of our current work. In contrast, based on Fig. 5, we 
expect the Square Kilometre Array (SKA) [80] to have sufficient precision for measurements of the 


non-Gaussian parameter. 
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FIG. 5. Energy-density fraction spectra h?Qew.o(v) (solid curves) and the variance of energy-density 
anisotropies h?[@(¢-+ 1)C;(v)/(2m)|!/? (dashed curves) versus the sensitivity of Square Kilometre Array 
(SKA) [80] (gray shaded area). 


V. CONCLUSIONS 


In this work, we studied the implications of the recent PTA datasets for the nature of primordial 
curvature perturbations and primordial black holes. In particular, we explored the influence of 
primordial scalar non-Gaussianity on the inferences of model parameters, and vice versa. By taking 
into account the impacts of primordial non-Gaussianity, we analyzed the current datasets and 
obtained the allowed parameter region for the primordial scalar spectral amplitude Ag ~ 107? — 1, 
the primordial non-Gaussian parameter —10 S Fynt < 10, and the mass of PBHs mpph ~ 1072 — 
1Mo. It is still important to stress the necessity of a full Bayesian analysis that will be conducted in 
the future. Even if the non-Gaussian parameter is considered, the PBH scenario was shown to be in 
tension with the NG15 data, except when we considered a sizable negative Fyi that can significantly 


suppress the abundance of PBHs. However, it is challenging to exclude the PBH scenario with 
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the current datasets, due to large uncertainties of formation models. Finally, we proposed that 
the anisotropies of SIGWs are a powerful probe for measurements of the non-Gaussian parameter 
Fyz, and conducted the first complete analysis of the angular power spectrum in the nano-Hertz 
band. In particular, we found that such a spectrum can effectively break some degeneracies of 
model parameters, particularly, the sign degeneracy of FNL. In addition, we explored detectablity 


of the anisotropies in SIGWs in an era of the SKA project. 
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Appendix A: Formulae for computing the SIGW energy density 


After the lengthy and tedious derivation in accordance with Refs. [31-33], the three terms in 


Eq. (4) can be exactly expressed as follows 


7 1 se 1 OO 1 
Ong) = al dt; [as J? (u1, 01, £ — 00) Taoa aa eM Aglaa) (Al) 
eo ON 
_ F2 2 o0 1 2 
20 (n, a) = St ff ati f ds; vs es thts Ohi 8) nA E a 
37 j= LJO —1 (uiv1U2v2) 
20 —_—_—_—_—_—_____ a raħie 
+ A dy12 cos 212 J (u1, v1, £ 4 00) J (u2, v2, £ —> 00) (A2) 
7 Aj (vq) A? (wi2q) E J: ak 
3 3 3 3 D 
U9 Wi2 ua ui 
5(2) FSL yrf [* i 
Qw (7, @) E 2472 I| at: f ds; vi 
A 
Qn? J? (u1, v1, £ > 00) 2 2 2 2 
7 A A A A 
(up vy ugVv2U303)2 g(U1U24) g(U1U29) g(U1U3q) g(U1u3q) 
2T r 
+ dyi2dyo3 cos 212 J (u1, 01, £ > 00) J (ua, v2, £ —> 00) (A3) 
0 


AZ (usq) A3(wisq) A2 (v3q) AZ (wo3q) A? (we3q) AZ(wiesq) 
4 3 3 3 3 t 3 3 , 
u3 W313 U3 w33 W953 W193 
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where we define £ = qn, Si = uj — Vi, ti = uj + vi — 1, and 


gig = SRP Heats + 2)(1 — 9B ty (ty + 2)C1 — 93) + FI silti + DII sft + 1)] (Ada) 


Wig = yV? + UF — Yij , (A4b) 


yo? + v3 + v? + yi2 — Y13 — Y23 - (A4c) 


W123 


The oscillation average of squared J(u, v, x — oo) has been given by Ref. [33] and the earlier works 


in Refs. [12, 13, 30, 31], i.e., 


— Tr.” 

J (tig, vi, £ > œ) J (uj, vj, £ —> 20) 

9 (1— s?) (1 — 8?) ti (ti + 2) ty (tj + 2) (8? + t? + 2ti — 5) (52 + 2 + 2t; — 5) 
E 8 (—si + ti + 1) 3 (s; + ti +1)’ (—s; + tj + 1) 3 (s; +t; +1)’ 


t? +2t; 2 
(2ta - s) m (FAEN ainat] 
7 
ete ee 
(Grd )m([EAT | ene) as 
2L 


+70 (ti — V3 +1) (tj — V3 +1) (s +8 + 2; -5) (434+ + 25-5) 


The formulae in this appendix can be used for numerically computing the energy density of SIGWs 


in a self-consistent way. 
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